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This  memorandum  describes  a  model  of  reactive  surveillance 
in  which  two  classes  of  contacts  occur.  A  single  search  vehicle 
attempts  to  maintain  some  localization  of  a  single  target  vehicle 
with  contacts  occurring  intermittently  as  a  Poisson  process  but 
only  when  the  target  is  within  a  certain  range  (exposure  disk)  of 
the  searcher.  An  external  surveillance  system  is  also  present; 
it  produces  contacts  as  a  Poisson  process  irrespective  of  the 
target  position.  Searcher  tactics  are  assumed  to  approximate 
optimal  search  based  on  an  assumed  circular  normal  target  dis¬ 
tribution  determined  solely  by  the  class  (searcher  or  external 
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The  remainder  of  this  memorandum  is  divided  into  sections. 

The  mathematical  model  and  possible  generalizations  are 
described  in  the  first  section.  Numerical  techniques  employed 
for  the  evaluation  of  limits  and  integrals  are  discussed  in  the 
second.  The  computational  algorithm  is  summarized  in  the  third 
section.  A  computer  program  listing  is  included  as  an  appendix. 
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Mathematical  Model 

Three  types  of  contacts  are  considered;  those  by  the 
searcher,  those  by  the  external  system  when  the  target  lies 
within  the  searcher's  exposure  disk,  and  those  by  the  external 
system  when  the  target  lies  outside  this  disk.  Probabilistically, 
the  sequence  of  contact  types  and  corresponding  contact  times  is 
described  as  a  Markov  renewal  process  in  which  the  distribution 
of  next  contact  type  and  recontact  time  depends  only  on  the  type 
of  most  recent  contact;  the  process  regenerates  with  each  contact. 
The  sequence  of  contact  types  is  a  Markov  chain.  The  intervals 
between  contacts  are  described  separately  as  inhomogeneous, 
continuous  time  parameter,  Markov  processes  with  different  time 
dependent  transition  functions  and  initial  probability  vectors 
according  to  the  type  of  most  recent  contact.  These  processes 
have  five  states  corresponding  to  the  three  types  of  contacts 
(absorbing  states)  and  target  locations  within  and  outside  of  the 
searcher's  exposure  disk  before  recontact.  The  assumed  circular 
normal  target  distribution  is  described  by  a  linear  datum  growth 
law  with  initial  datum  area  determined  by  contact  type  and  datum 
growth  rate  determined  by  target  motion  statistics. 

These  concepts  are  intended  to  model  reactive  surveillance 
by  ANVs;  the  external  contacts  might  correspond  to  those  generated 
by  SOSUS  or  by  all  other  surveillance  systems  combined.  However, 
the  concepts  generalize  to  situations  involving  differing  numbers 
of  contact  types,  transition  functions,  and  datum  growth  laws. 

The  principal  measure  of  effectiveness  derived  in  this  memo¬ 
randum  is  distribution  of  actual  target  range  from  the  expected 
(by  the  searcher)  target  location  at  a  random  time  instant.  This 
has  obvious  application  to  anti-SLBM  defense.  Many  other  such 
measures  could  be  derived  from  the  information  generated  by  the 
model. 
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The  five  states  in  the  interval  following  a  contact  are 
numbered  and  described  as  follows: 


State  1  -  Recontact  by  the  external  surveillance 
system  has  occurred  while  the  target  is 
exposed  to  contact  by  the  searcher. 

State  2  -  Recontact  by  the  external  system  has 

occurred  while  the  target  is  unexposed 
to  the  searcher. 

State  3  -  Recontact  by  the  searcher  has  occurred. 

State  4  -  Recontact  has  not  yet  occurred  and  the 
target  is  exposed. 

State  5  -  Recontact  has  not  yet  occurred  and  the 
target  is  unexposed. 


The  three  types  of  contacts  have  the  same  numbering  as  the  first 
three  states. 


Let 


1,2, 3, 4, 5 
1,2,3 


and 


Pj(t)  -  (Plj(t),-**,P5j(t)),  (j  -  1,2,3) 

where 


t  “  time  since  last  contact. 


Pij(t) 


■  Probability 


ity  ^ 


process  in 
state  i  last 
at  time  t  contact  of 
after  last  type  j  , 
contact 


Three  separate  Markov  processes  are  obtained  for  j  ■  1,2,3.  By 
definition, 


P^o)  -  (0,0,0, 1,0) 


and 


P2(0)  -  (0,0, 0,0,1). 


Since  the  target  must  be  exposed  for  contact  by  the  searcher  to 
occur , 


P3(0)  -  (0,0, 0,1,0). 


The  model  uses  three  absorbing  states  corresponding  to 
contact  types  although  the  use  of  only  two  (external  and  by 
searcher)  would  appear  more  natural.  This  is  done  so  that  the 
initial  probability  vectors,  P-^(O),  can  be  specified  a  priori. 

In  adapting  the  model  to  a  more  general  situation  in  which  there 
are  N  (in  this  case  two)  classes  of  contacts  and  M  (in  this  case 
also  two)  intervening  states,  as  many  as  MN  (in  this  case  three 
sufficed)  absorbing  states  might  be  required. 


Let 


5(A) 


unconditional  probability 
that  target  is  exposed  to 
searcher  during  search  of 
datum  area  A 


and 


*(A) 


rate  at  which  unexposed  targets 
enter  exposure  disk  during 
search  of  datum  area  A. 


d 


dT 


Probability 


target  exposed 
at  time  t  +  t 


target 
unexposed 
at  time  t 


1 


T-0. 


5 


"Datum  area"  is  defined  by 

A  -  6tto2  (1) 

2 

where  o  is  the  parameter  characterizing  the  assumed  circular 
normal  target  distribution.  This  is  the  area  of  a  bounded  region 
which,  if  it  contained  the  target  with  a  uniform  location  distri¬ 
bution,  would  require  the  same  expected  effort  to  detect  as  is 

2 

required  by  a  circular  normal  distribution  with  parameter  o 
(see  reference  Ca3).  While  in  actuality  the  quantities  §  and  » 
depend  on  geographic  details  and  are  not  functions  of  A  alone, 
the  following  formulas  are  reasonable  approximations: 

5(A)  -  1  -  exp(-n  R£2/a)  (2) 

and 


i(A)  -  2ReVs/A, 

where 


(3) 


R£  ■  radius  of  exposure  disk,  and 
V  -  searcher  patrol  speed. 

Formula  (2)  subsumes  searcher  allocation  of  effort  and  non¬ 
regularity  in  the  search  region.  Formula  (3)  is  familiar  and 
2 

standard  for  ttR  «  a,  it  is  shown  below  that  it  is  sensible  for 
small  A  as  well. 

Consider  an  arbitrary  Re  and  A.  If  the  target  is  not 
exposed  to  detection,  it  is  confined  to  a  region  of  expected  area 
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Let  “2  be  the  speed  at  which  the  boundary  of  the  disc  of  exposure 

sweeps  through  this  unexposed  area  and  let  W  be  effective  sweep 
2 

width.  For  fTRg  «  A,  the  entire  boundary  intersects  A, 

*S  ■  V  and  W  *■*  2Rg.  For  TTRg^  »  A,  only  a  small  portion  of  the 
boundary  intersects  A,  and  speed  perpendicular  to  that  segment 
is  V_( 2/n).  But  with  probability  .5,  that  segment  is  sweeping 
into  exposed  area*  i.e.,  no  new  area  is  exposed.  Thus  S  **  S/tt 
and  W  the  length  of  boundary  intersecting  A. 

Now  if  the  target  is  located  in  the  unexposed  area*  the 
rate  of  entry  into  the  disc  of  exposure  is 

SW 

A  exp  (-ttRe  /A) 

Equating  this  to  (3),  one  has 


which,  for  ttRe  «  A,  gives  the  expected  result,  W m  2Rg.  As  A 
decreases  in  size*  the  right-hand  member  of  (4)  decreases*  and 
can  be  bounded  by 


W*  (t5t)(£)  (f^2  exp  (-.5))  i 


Taking  *2  -  S/n *  a  presumed  minimum  value, 


(4) 
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(2.694), 


W  * 


a  very  reasonable  result.  Thus  formula  (3)  appears  sensible 
for  small  A,  despite  the  inverse  power  dependence. 

The  assumed  linear  datum  growth  law  after  contacts  of  type 

j  is 


Aj(t)  "  Aj(0)  +  pt  (5) 

where 

Ai(t)  ■  the  datum  area  at  a  time  t  after  a  contact 
J  of  type  j,  and 

p  ■  datum  growth  rate, 

for  j  “  1,2,3. 

Since  type  1  and  type  2  contacts  both  come  from  the  same  external 
surveillance  system,  the  resulting  initial  datum  is  the  same  for 
each  (A^O)  -  A2(0)).  The  growth  rate  p  is  determined  by  target 
motion  statistics.  If  the  target  is  moving  in  an  unbiased  random 
walk  with 


-  target  speed,  and 

C  “  mean  time  between  target  course  changes, 


it  may  be  shown  that  the  target  location  distribution  is  approxi¬ 
mately  circularly  normal  with  density  in  polar  coordinates  fr,0). 


f(r,e) 


2TTOZ(t) 


exp 


(6) 
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¥**-■ 


where 

o 2(t )  -  Vt2  C2  [t/C-  (1-  exp(-t/C))3. 
Thus  for  large  t, 

So2(t)  9 


9 

and,  since  A  *  6tto  , 

p  -  6TTVt2C  .  (7) 

Formula  (7)  is  used  to  compute  p  in  applications. 


Let 


a  -  external  recontact  rate 


d 


&T 


|  Probability 


external 
recontact 
time  t  +  t 


by 


no  recontact 
by  time  t 


» 

T-0 


X  ■  searcher  recontact  rate  for  an  exposed  target 


d 

dr 


Probability 


searcher 
recontact 
t ime  t  +  t 


1# 


by 


target  exposed 
at  time  t  and 
no  recontact 
by  time  t 


» 


T-0 
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0i(t)  ■  rate  at  which  unexposed  targets  enter  the 
J  exposure  disk  after  a  type  j  contact 


Probability 


target 
exposed  at 
time  t+t 


target  unexposed] 
at  time  t  and 
last  contact  of 
type  j 


* 


T-0 


and 


Vj(t) 


rate  at  which  exposed  targets  leave  exposure  disk 


<  Probability 


target 
unexposed 
t ime  t  +  t 


at 


target  exposed 
at  time  t  and 
last  contact  of 
type  j 


T=0 


where 


t  ”  time  since  last  contact, 


for  j  ®  1,2,3. 


It  is  assumed  that  the  rates  a  and  X  are  actually  constant 
and  that  the  rates  pj  and  y  j  are  functions  of  t  alone.  The  model 
requires  these  assumptions,  which  seem  quite  reasonable.  It 
follows  that 


dpj(t) 

dt 


PJ(t) 


0 

0 

0 

a 


0 

0 

0 

0 


0 


a 


0  0  0 

0  0  0 

0  0  0 

X  -ta+\+Yj(t)]  Y  j  (t) 

0  0j(t)  “[a>+p  j(t)] 


(8) 
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Now 


Pj(t)  -  ’(Aj (t ) ) 


_  2REVs 
Aj(t) 


W 


by  formula  (3).  For  the  case  a  -  \  «  0,  P^j(t)  "  §(Aj(t))  and 
the  fourth  equation  in  (8)  becomes 


? ' (Aj(t) )Aj ' (t) 


-Yj(t)S(Aj(t))+Pj(t)(l-5(Aj(t))). 


Given  a  form  for  5  this  can  be  solved  for  Yj  in  terms  of  pj  and 
Aj.  Using  formula  (2)  for  §, 


Yj(t) 


'  exp(-HRE2/A1(t)) 
l-exp(-T'RE2Aj(t)) 


9j(t) 

- 


+ 


^ Aj (t)  “J 
at 


CAJ(t)32 


exp(-TTRE2/A.(t)) 

1  -  exp(-TTRE  Aj(t)) 


Pj(t)  + 


CAj(t)32 


(10) 


System  (8)  could  be  easily  modified  to  discard  external 
contacts  occurring  at  times  when  their  information  is  inferior 
to  that  already  available  from  the  last  contact.  This  would 
entail  replacing  the  term  a  in  (8)  by 

0  if  A.(t)  <  A1(0) 
a  otherwise 

This  modification  would  have  little  effect  in  most  cases. 
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Let 


BtJ  -  Pjjtt)  for  i,J  -  1,2,3. 


It  may  be  shown  that  since  a  >  0,  the  matrix  tBjj3  is  well 
defined*  has  all  positive  entries,  and  all  column  sums  one. 

[Bij]  is  the  transition  matrix  of  the  Markov  chain  characterizing 
the  sequence  of  contact  types.  It  has  a  unique  eigenvector 
summing  to  one  with  eigenvalue  one;  this  is  the  vector  of  steady 
state  probabilities.  Let 


<91>92’93>  ’ 


steady  state  probability  vector  for 
Markov  chain  characterizing  sequence 
of  contact  types 


(ID 


If 


®13^®22  ""  ^)  “ 


23  12 


(Bn  -  1)(B22  -  1)  -  B12B21 


Define 


®23^®11  “  ^)  “ 


13  21 


(Bj^  ”  1)(B^^  -  1)  -  B-|  <jBr 


and 


'22 


12  21 


*3  "  ‘I 


then 


(12) 


e,  - 


± 


j  *l+i2+i3 


for  j  -  1,2,3. 


Qj(t) 


Probability 


recontact  within 


recont 

time  t 


last  contact 
of  type  j  J 
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and 


A  j(t) 


Qj(s)]  ds  . 


Then 


Qj(t)  "  i«!  Pij(t)  * 


(13) 


(14) 


Integrating  by  parts, 


Conditional  expectation 
of  duration  of  portion 
of  interval  between 
contacts  where  time 
since  last  contact  <  t 
given  that  last  contact 
„of  type  j 


I’*  (  s  if  s  s  t' 

J0  (t  if  s^t, 


-  fij(t). 


(15) 


In  particular,  defining 


lim 
t  «*  • 


n:<c)- 


(16) 


Lj 


Expected  value  of 
time  until  recontact 
after  contact  of 
type  j 


(17) 


This  limit  certainly  exists  since  a  >  0  and  (8)  imply  that 
l-Qj(t)  is  bounded  above  by  exp(^xt). 
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for  which  the  initial  state  is  assumed  known  .  Let 

T  .  a  random  time  instant,  uniformly 

[distributed  on  [0,T^]  J  * 

i  m  the  type  of  the  last  contact  preceding 
J1  T  (0  if  there  is  no  such  contact)  J  * 

*  m  the  type  of  the  first  contact  following"]  . 

J2  T  (0  if  there  is  no  such  contact)  J  *  ana 

rp  „  time  elapsed  at  T  since  last  contact 
le  (0  if  there  is  no  such  contact)  * 

The  random  variables  J2*  and  Tg  all  have  limiting  dis¬ 
tributions  as  which  are  independent  of  the  initial  state. 

These  distributions  are: 

Probability  tjx  -  jJ  -  ^  (by^ll'.^n)!  and  (20))’ 


Probability  [j2  -  j]  =  iii  (^({ij, 


and  (20); 


(23) 


Probability  [Te  <  t] 


S  8knk(t) 


r  \by  (11),  (15),  and  (20)/’ 


Now  define 


F(V)  -  Probability  ?]  ,  and 


Aj'^Y)  - 


0  if  Y  s  Aj (0) 

unique  t  such  that  Aj(t)  -  Y  otherwise 


Here  the  model  is  regarded  as  a  five  state  process  with  continu¬ 
ous  time  parameter  in  which  the  three  states  corresponding  to 
contacts  are  attained  only  for  discrete  instants. 


15  - 


Then  by  (24), 


2  0^k(Ak’1(Y)) 
F(Y)  -  — - 

r 


With  an  assumed  circular  normal  datum  with  parameter  o  ,  formula 
(6)  implies 


target  within 
Probability  distance  R  of 
datum  center 


1  -  exp(-R2/2o2) 


Defining 


Unconditional  steady  state  probability 
G(R)  *  that  target  is  within  distance  R  of  , 
datum  center 


G(R) 


Jf*  0  dF(6Trx) 

I  [  1  -  exp(-Rz/2x)3  -  dx  ,  using  (1), 

0  dx 


Integrating  by  parts  and  changing  variables, 


G(R)  -  3ttR 


2  /**  exp(-3TTR2/U)F(U) 


Using  formulas  (25)  and  (5)  for  the  functions  F  and  Aj  and 
changing  variables  again  leads  to 


G(R) 


£  i  Lf 

r  \^i\  kJ0 


expC-3TTR2/(A^(0)  +pt)3 
~  (A^O+pt)* 


nk(t)dt[.  (26) 
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Evaluation  of  Limits  and  Integrals 


System  (8)  is  evaluated  numerically  by  the  Runge-Kutta 
method  for  j  ■  1,2,3.  The  vectors  P^(t)  -  (P^ j (t) , *  *  * , j(t)) 
and  dP-Vdt  -  ' (t) ,  *  *  *  ,Pg  j ' (t))  are  available  explicitly  at 

each  increment.  The  evaluation  is  carried  out  over  the  interval 
C0*  Tmax]»  where  Tmax  is  chosen  so  that  the  recontact  probability, 
QjCTjj^), will  be  quite  large  (say*  .95)  for  each  j.  The  contact 
type  transition  probabilities  are 


Bij  - 11-- 


Since 

3 

S  B..  -  1 
i-1  1J 


necessarily,  a  convenient  and  reasonable  approximation  is 

Bu~  VW*  a1J’(Tmix)  (x- J,pu<w)  •  (27) 

The  functions  0 j C t )  were  defined  as 

flj(t)  “  f  [1-Qj(s)3ds 
•'0 
3 

where  Qj(s)  -  Pjj(s).  Since  values  of  the  Pjj  and  their 
derivatives  are  explicitly  available  at  fine  increments  from  the 
Runge-Kutta  evaluation  of  system  (8),  the  0^  can  be  obtained  at 
the  same  time  to  reasonable  accuracy  by  trapezoidal  integration 
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1 


for  times  in  the  interval  (0,  .  The  limits 

Lj  -  V0 

are  also  required.  Approximations  to  the  Lj  are  obtained  by 
assuming  that  the  probability  of  no  recontact  becomes  a  negative 
exponential  after  a  long  time.  This  assumption  is  reasonable, 
particularly  if  datum  growth  really  ceases  after  a  long  time. 
More  precisely  it  is  assumed  that  for  t  2  Tfnax, 


1-Q,(t) 

- ] -  *  exp[-*.(t  -  T  )3 

1_  (T  \  J  tnax 

max/ 


for  some  constant  i y  This  leads  to  the  approximation 


x  ^-Ql^max^2 

L,  ~  0  ,(Tmov)  + - 3 — — -  , 

J  j  max  0  '  fT  ) 


j  '  max' 


which  is  calculable  since  fij,  Qj,  and  Q ^ ’  are  all  available  for 


t  ~  Tmax* 


The  integrals  in 


G(R) 


3TTPR" 


I T  k-1 


tM 


exp[-3TTR2/(Ak(0)  +  p  (t)3 

■2 - «k(t)dt 


(Ak(0)  +pt)' 


(28) 


must  be  evaluated  for  a  number  of  values  of  R.  The  portion  of 
each  integral  in  the  interval  (0,  T_ov.]  is  evaluated  from  a  table 
of  values  of  (t,fik(t))  which  was  saved  from  the  Runge-Kutta 
integration.  The  increment  size  changes  several  times  in  the 
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table*  Simpson's  rule  is  used  separately  on  each  Interval  with 
a  constant  increment  size;  trapezoidal  integration  is  used 
wherever  an  odd  subinterval  is  left  over.  These  integrals  have 
substantial  "tails"  in  the  interval  [T  •].  These  are 

IUclX. 

estimated  by  assuming  that  for  t  >  Tmax, 

nk(t)  ~  Lk,  and 
Ak(0)  +  pt  -  pt. 


This  leads  to  the  approximation 


3ttpIT  3 

G(R)  *  -  S  l  0, 

I T  k-1 


l 


rmax  exp[-37TR2/(Ak(0) +pt)] 


(Ak(0)+pt) 


2 - H^Ct)  dt 


+  l-expC-STTR^/pT^)  . 


(29) 
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Summary  of  Computational  Algorithm 

The  algorithm  is  embodied  in  a  computer  program  which  is 
included  as  Appendix  A.  Inputs  are: 

V.  ■  searcher  patrol  speed  (nm./hr.), 

R£  “  radius  of  exposure  disk  (nm.), 

a  -  external  recontact  rate  (hr.”^)  when  target 
is  exposed, 

X  ■  searcher  recontact  rate  (hr.”*0, 

o 

p  *  datum  growth  rate  (nm./hr.), 

Ai(0)  ■  initial  datum  size  after  external  contacts 
(nm.  ),  and 

A3(0)  -  initial  datum  size  (nm.  )  after  searcher 
recontacts. 


The  systems  of  linear  first-order  differential  equations 
(8)  are  Integrated  over  the  interval  [0,  Tmny]  by  a  library 
Runge-Kutta  routine  for  j  ■  1,2,3.  The  rates  Pj(t)  and  Yj(t) 
are  calculated  by  formulas  (9)  and  (10)  using  the  datum  growth 
law  (5).  Values  of  Qj(t)  and  fij(t)  are  calculated  simultaneously 
as  discussed  earlier  using  formulas  (13)  and  (14).  Values  of  t 
(time  since  last  contact)  and  Oj(t)  are  retained  for  later  use 
for  a  number  of  increments.  The  reason  for  retaining  the  time 
values  is  that  the  Runge-Kutta  routine  determines  and  adjusts  its 
own  step  size,  defying  user  control.  Values  of  t,  Pjj(t), 
p2j(t)»  P3j(t),  P5j(t),  Qj(t),  and  Oj(t)  are  printed  for 

a  number  of  values  of  t. 


ft 


The  limits  and  Lj  are  calculated  by  formulas  (27)  and 
(28).  The  steady  state  contact  type  probabilities  are 
calculated  by  formula  (12).  The  Bjj,  Lj,  and  6j  are  then 
printed  out. 

Values  of  the  range  distribution  G(R)  are  calculated  from 
the  tables  of  t,  fij(t)  an<*  printed  for  the  values  R  **  5,10, 

15, ’**,200  nm. ,  as  discussed  earlier  using  formula  (29). 


Peter  S.  Shoenfela 


Reference  [a]:  L.  D.  Stone,  Theory  of  Optimal  Search.  Academic 
Press,  New  York,  1975. 
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APPENDIX  A 


COMPUTER  PROGRAM  LISTING 


/ 


C 

5 

10 

20 

C 

C 


30 


50 


100 

110 

C 

C 


200 

C 


250 


DIMENSION  PRMT  <  5 )  f DERP  <5)fP<5)f  AUX  <  8  f  5 ) 

DIMENSION  PH<3f3>  fDPD<3f3> fT<3f400) f0MEGA<3f400> 

DIMENSION  OMINF < 3 ) f THETA (3)fB(3f3)fN(3) 

REAL  LAMDDA 

COMMON  VS » ALPHA f R r LAMBDA r  RHO f AO  1 1 A02 f J » PD » DPD r T V OMEGA t 
2  OMINF fNfTLAST fPDLAST 
DATA  PRMT/O. r 400. f 0.0625 f 0.0001 fO./ 

DATA  NDIM/5/ fTMAX/300 • / 

READ  PARAMETERS 
TYPE  10 

FORMAT ( 1 X  f  ' VS  t RE t ALPHA  f  LAMBDA  f  RHO  f  AZER01 f  AZER03? ' / > 

ACCEPT  20  f  VS  f  R  f  ALPHA  f  LAMBDA  f  RHO  f  A01 f  A02 
FORMAT <7G> 

EXTERNAL  FCTfOUTP 

INTEGRATE  O.D.E.  SYSTEM  FOR  EACH  CONTACT  TYPE 
BY  RUNGE  KUTTA  METHOD  USING  SUBROUTINE  RKGS 
PRMT  <  2 ) =TMAX 
DO  100  J=1f3 
PRMT  <5)=0. 

DO  30  N=1 f  5 
P ( K ) =0 . 

DERP<K>=0.2 

IF< J.NE.2)P(4>-1. 

IF<J.EQ.2)P<5>=1. 

N( J)=0 
OMINF  <  J ) =0 • 

TLAST =0 . 

PDLAST =0 • 

TYPE  50  f  J 

F0RMAT<20Xf' TYPE ' f 13  f  '  CONTACT ' //4X f 'T' f  SXf'PI'f  4Xf 
2  '  P2 '  f  4X  f  ' P3 ' f  4 X f '  F"  4 ' r  4Xf'P5'f  4X f ' PD ' f 4X f ' OMEGA '// > 

CALL  RKGS ( PRMT  f F* f  DERP  f  ND I M f  IMLF f FCT  f OUTP f  AUX ) 

TYPE  IIOfIHLF 

FORMAT ( IXf '  IHl.F='  f  I4//> 

FIND  LIMITS  TO  OBTAIN  TRANSITION  MATRIX  B(IfJ)  AND 
MEAN  HOLDING  TIMES  OMINF<J) 

DO  200  J=1 , 3 

PDTOT =PD ( J  f 1 ) +PD <  J  f  2 ) +PD  <  J  f  3 ) 

DPDTOT =DF‘D  ( J  f  1 )  +DPD  (  J  f  2 )  +DPD  <  J  f  3  > 

OMINF < J)=OMINF< J)+( 1. -PDTOT >**2/DPDT0T 
DO  200  1  =  1 f  3 

B  < I f  J ) =PD  <  J  f I ) +DPD ( J  *  I ) * ( 1 . -PDTOT  > /DPDTOT 

FIND  EIGENVECTOR  OF  B(IfJ)  TO  GET  STEADY  STATE  PROBABILITIES 
D=(B<1f1)-1. >*<B(2f2>~1. )-B<1f2>*B(2f1) 

THETA < 1)»(B<1f3)#<B(2f2)-1 . > -B < 2 f 3 > *B < 1 f 2 >  > /D 
THETA  <  2  )  =  <B<2f3)JK(B(1f1>-1.)-E«(1f3)#B<2f1)>/D 
THETA<3)=-1 ♦ 

SUM=THETA<  1  )+THETA(  2 )  +THETA<  3  ) 

DO  2-0  J=1f3  r BRA-1 

THETA  ( J>=THETA<  J)/SUM  IHIS  ^  lOUDC  _ - ^ 


C  PRINT  THESE  RESULTS 

TYPE  260f0MINFf THETA*  <  <B(If J) » J-l»3>  *1=1*3) 

260  FORMAT ( //20X *  ' MEiiAN  HOLDING  TIMES'// 

2  3<3X*Er13.5> //20X  *  '  STEADY  STATE  PROBABILITIES V/3<3X»F13f5>// 

3  20X * ' TRANSITION  MATRIX ' //3< 3< 3X *F13 .5) / > ) 

C  OAl.CtM  ATI:  RANOL  HI  SIR  I  HUT  ION  AND  PRINT 

TYI  L  2/0 

2/0  FORMA  I  <  // 1 X *  '  RANGE '  *  SX  r  '  CUM .  PROD .  '  // ) 

DO  600  K«J1  *  40 
RR»3»*K 
PROD  -*0. 

D0  300  J«l»3 
SUM-0. 

NN=0 

T0LD=0. 

290  IF<NN.EQ.N( J) )G0  TO  400 

H1=T ( J  *  NN+1 ) -TOLD 

IF<NN+2.LE.N< J>  >H2*T< J*NN+2>-T< J*NN+1> 

IF < M2 . NE . HI . OR . NN+2 ♦ GT . N < J) )G0  TO  300 

SUM=SUM+H1*<ELT  (NN?RR)+4» KELT  <NN+1 *RR  >  +ELT  <NN+2*RR>  >/3. 

TOLD=T  <  J  *  NN+2 ) 

NN=NN+2 
GO  TO  280 

300  SUM=SUM+H1*<ELT<NNfRR)+ELT<NN+1fRR>  >/2. 

TOLD-T <  J  * NN+1 ) 

NN=NN+1 
GO  TO  280 

400  PR0B=PR0B+THETA< J)*SUM 

500  CONTINUE 

PROB-9 . 424778*RH0*RR*RR*PR0B/  <  THETA  <  1  >  JKOMINF  <  1  >  + 

2  THETA<2)*OMINF(2)+THETA(3>#OMINF(3) )+ 

3  1 . -EXP ( -9 . 424778*RR*RR/ < RHO*TMAX ) ) 

TYPE  510 »RR» PROD 

510  FORMAT  <1XfF6.1*2XfF7.3) 

600  CONTINUE 
GO  TO  5 
END 

FUNCTION  ELT < I »RR> 

C  EVALUATES  INTEGRAND  FOR  RANGE  DISTRIBUTION  CALCULATION 

DIMENSION  PD  <  3  *  3  >  f  DPD  <  3 * 3 ) *  T ( 3  * 400 ) t OMEGA ( 3 » 400 >  » 

2  OMINF ( 3 ) f THETA <3>fB(3f3)fN(3) 

REAL  LAMBDA 

COMMON  VS  f  ALPHA  f  R f  LAMBDA  f  RHO f  AO 1 f  A02 f  J  f PD  f  DPD f  T  f  OMEGA  f 
2  OMINFfNfTLASTfPDLAST 
A0=A01 

IF( J.EQ.3)A0=A02 

IF<I.GT.0» AND «I.LE.N(J))GO  TO  20 
ELT  ==0  • 

GO  TO  30 
20  TT=T(JfI) 

ELT=EXP<-9.424778*RR*RR/<A0+RH0*TT>  >*OMEGA( JfI)/ 

2  <A0+RH0*TT>**2 
30  RETURN 

END 


oi  r  j 


I - 

i 


C 


C 


10 


0 

0 


70 


* 


SUBROUTINE  FCKTT. P.DERP) 

EVALUATES  DERIVATIVES  FOR  RUNGE-NUTTA  ROUTINE 
DIMENSION  P<5>,DERP<5> 

DIMENSION  PD<  3 » 3 ) »  DPD ( 3 » 3 )  t  T  ( 3 » 400  >  t OMEGA <  3 1 400 >  t 
2  OMINF <  3 > » THETA ( 3)»B(3»3)»N<3) 

REAL  LAMBDA 

COMMON  VS » ALPHA  >  F< » LAMBDA ,  RHO » AO  1 » A02 »  J ,  PD  *  DPD  r  T » OMEGA » 
2  OMINF »N»TLAST »PDLAST 
A(X)=AO+RHO*X 
BETA<  X  )=2 . #R#VS/A ( X ) 


PP(X)=EXP<-3.141593*R*R/A<X> ) 

GAMMA<X)=<PP(X)/(1 .-PP<X> ) >*<BETA<X)+3.141593*R*R*RH0/A<X>**2> 
A0=A01 

IF (J.EQ.3)A0=A02 
DERP<1)=ALPHA*P<4> 

DERP  <  2 )  ==ALPHA*P  <  5 ) 

DERP  <  3 )  =LAMBDA*P  ( A  ) 

DERP  <  4  )  =-  ( ALPHA+LAMBDA-f  GAMMA  ( TT  >  >  *P  ( 4  )  +  BETA  ( TT )  *P  ( 5 ) 

DERP  <  5 ) =GAMMA <  TT ) *P ( 4 ) - ( ALPHA+BETA ( TT ) ) *P < 5 > 

RETURN 

END 


SUBROUTINE  OUTP< TT » P , DERP » IHLF » NDIM * PRMT ) 

OUTPUT  ROUTINE  FOR  RUNGE-KUTTA  INTEGRATION 
REAL  LAMBDA 

DIMENSION  P  ( 5 ) » DERP ( 5 )  » PRMT < 5 ) 

.  DIMENSION  PD<3»3) »DPD(3»3)  »T (3»400) »QMEGA(3r400>  t 
2  OMINF (3) >THETA(3) »B(3»3) »N<3> 

COMMON  VS  >  ALPHA » R » LAMBDA » RHO  t  A01  r  A02 » J»  PD » DPD  r  T  r  OMEGA r 
2  OMINF »N»TLAST rPDLAST 
,  DO  10  1=1 , 3 
F'D <  J*  I  )=P<  I ) 

DPD  <  J » I ) =DERP  < I ) 

PDD=P  < 1 ) +P  <  2  > +P  <  3 ) 

OMINF  <  J ) =OMINF  <  J )  +  < TT-TLAST ) *< 1 . -(PDD+PDLAST) /2 .  > 

TLAST =TT 
PDLAST=PDD 

IF<AMOD<TT, ,25). GE. 0.0001 )G0  TO  50 
IF < AMOD ( TT » 2 . ) . GE ♦ 0 . 0001 . AND . TT , GT . 4 . >G0  TO  50 
IF<AM0D<TT»8.).GE. 0.0001. AND. TT.GT. 20. >G0  TO  50 
IF ( AMOD < TT f 16 . ) • GE. 0 • 0001 • AND . TT . GT , 100 , ) GO  TO  50 
TYPE  20 f  TT r P » PDD f OMINF < J > 

FORMAT  <lXrF7.2»6<lX»F5.3)flX»F6.2) 

IF(AMOD(TTf ,0625) .GE. 0.0001 )G0  TO  70 

IF(AMOD<TTr 1. ) .GE. 0.0001. AND. TT.GT. 4. >G0  TO  70 

N ( J )=N  <  J ) +1 

T(J?N(J>) =TT 

OMEGA ( J »N <  J > ) ®OMINF ( J ) 

IF (PDD.GT . 1 . )PRMT (5)  =  1 . 

RETURN 

END  lr 
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in  um.vi;:  a  system  ok  first  order  ordinary  differential 
i;i.iiwiiimi  with  o i vi  in  inihae  values. 

USAGE 

CALL  RKOS  (PRMT » Y .DERY .NDIM. IHLF.FCT .OUTP.AUX) 
PARAMETERS  FCT  AND  OUTP  REQUIRE  AN  EXTERNAL  STATEMENT. 


DESCRIPTION  OF  PARAMETERS 

PRMT  -  AN  INPUT  AND  OUTPUT  VECTOR  WITH  DIMENSION  GREATER 


OR  EQUAL  TO  5 »  WHICH  SPECIFIES  THE  PARAMETERS  OF 
THE  INTERVAL  AND  OF  ACCURACY  AND  WHICH  SERVES  FOR 
COMMUNICATION  BETWEEN  OUTPUT  SUBROUTINE  (FURNISHED 


BY  THE  USER)  AND  SUBROUTINE  RKGS •  EXCEPT  PRMT(5) 

THE  COMPONENTS  ARE  NOT  DESTROYED  BY  SUBROUTINE 
RKGS  AND  THEY  ARE 

PRMT ( 1 ) -  LOWER  BOUND  OF  THE  INTERVAL  ( INPUT). 

PRMT ( 2 ) -  UPPER  BOUND  OF  THE  INTERVAL  (INPUT). 

PRMT ( 3 ) -  INITIAL  INCREMENT  OF  THE  INDEPENDENT  VARIABLE 
(INPUT)  . 

PRMT ( A ) -  UPPER  ERROR  BOUND  (INPUT).  IF  ABSOLUTE  ERROR  IS 
GREATER  THAN  PRMT (A).  INCREMENT  GETS  HALVED. 

IF  INCREMENT  IS  LESS  THAN  PRMT<3>  AND  ABSOLUTE 
ERROR  LESS  THAN  PRMT ( A ) /50 .  INCREMENT  GETS  DOUBLED. 
THE  USER  MAY  CHANGE  PRMT ( A  >  BY  MEANS  OF  HIS 
OUTPUT  SUBROUTINE. 


PRMT ( 5 ) - 


Y 


DERY 


NDIM 

IHLF 


NO  INPUT  PARAMETER.  SUBROUTINE  RKGS  INITIALIZES 
PRMT (3) =0.  IF  THE  USER  WANTS  TO  TERMINATE 
SUBROUTINE  RKGS  AT  ANY  OUTPUT  POINT.  HE  HAS  TO 
CHANGE  PRMT ( 5 )  TO  NON-ZERO  BY  MEANS  OF  SUBROUTINE 
OUTP.  FURTHER  COMPONENTS  OF  VECTOR  PRMT  ARE 
FEASIBLE  IF  ITS  DIMENSION  IS  DEFINED  GREATER 
THAN  5.  HOWEVER  SUBROUTINE  RKGS  DOES  NOT  REQUIRE 
AND  CHANGE  THEM.  NEVERTHELESS  THEY  MAY  BE  USEFUL 
FOR  HANDING  RESULT  VALUES  TO  THE  MAIN  PROGRAM 
(CALLING  RKGS)  WHICH  ARE  OBTAINED  BY  SPECIAL 
MANIPULATIONS  WITH  OUTPUT  DATA  IN  SUBROUTINE  OUTP. 
INPUT  VECTOR  OF  INITIAL  VALUES.  (DESTROYED) 

LATERON  Y  IS  THE  RESULTING  VECTOR  OF  DEPENDENT 
VARIABLES  COMPUTED  AT  INTERMEDIATE  POINTS  X. 

INPUT  VECTOR  OF  ERROR  WEIGHTS.  (DESTROYED) 

THE  SUM  OF  ITS  COMPONENTS  MUST  BE  EQUAL  TO  1. 
LATERON  DERY  IS  THE  VECTOR  OF  DERIVATIVES.  WHICH 
BELONG  TO  FUNCTION  VALUES  Y  AT  A  POINT  X. 

AN  .INPUT  VALUE.  WHICH  SPECIFIES  THE  NUMBER  OF 
EQUATIONS  IN  THE  SYSTEM, 

AN  OUTPUT  VALUE.  WHICH  SPECIFIES  THE  NUMBER  OF 
BISECTIONS  OF  THE  INITIAL  INCREMENT.  IF  IHLF  GETS 
GREATER  THAN  10.  SUBROUTINE  RKGS  RETURNS  WITH 
ERROR  MESSAGE  IHLF^U  INTO  MAIN  PROGRAM.  ERROR 
MESSAGE  IHLF- 12  OR  IHLF-13  APPEARS  IN  CASE 
PRMT ( 3 ) -0  OR  IN  CASE  SIGN ( PRMT ( 3 )). NE . SIGN (PRMT<2>- 


PRMT ( 1 ) )  RESPECTIVELY. 
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FCT  -  THE  NAME  OF  AN  EXTERNAL  SUBROUTINE  USED.  THIS 

SUBROUTINE  COMPUTES  THE  RIGHT  HAND  SIDES  DERY  OF 
THE  SYSTEM  TO  GIVEN  VALUES  X  AND  Y.  ITS  PARAMETER 
LIST  MUST  BE  X  *  Y*  DERY .  SUBROUTINE  FCT  SHOULD 
NOT  DESTROY  X  AND  Y. 

OUTP  -  THE  NAME  OF  AN  EXTERNAL  OUTPUT  SUBROUTINE  USED. 

ITS  PARAMETER  LIST  MUST  BE  X»Y *DERY * IHLF »NDIM»PRMT • 
NONE  OF  THESE  PARAMETERS  (EXCEPT *  IF  NECESSARY » 

PRMT  ( 4  >  * PRMT ( 5 )  *  .  .  .  )  SHOULD  BE  CHANGED  BY 
SUBROUTINE  OUTP.  IF  PRMT<5>  IS  CHANGED  TO  NON-ZERO» 
SUBROUTINE  RKGS  IS  TERMINATED. 

AUX  -  AN  AUXILIARY  STORAGE  ARRAY  WITH  8  ROWS  AND  NDIM 
COLUMNS. 

REMARKS 

THE  PROCEDURE  TERMINATES  AND  RETURNS  TO  CALLING  PROGRAM*  IF 

(1)  MORE  THAN  10  BISECTIONS  OF  THE  INITIAL  INCREMENT  ARE 
NECESSARY  TO  GET  SATISFACTORY  ACCURACY  (ERROR  MESSAGE 
IHLF*=1 1  >  * 

(2)  INITIAL  INCREMENT  IS  EQUAL  TO  0  OR  HAS  WRONG  SIGN 
(ERROR  MESSAGES  IHLF=12  OR  IHLF=13>» 

(3)  THE  WHOLE  INTEGRATION  INTERVAL  IS  WORKED  THROUGH* 

(A)  SUBROUTINE  OUTP  HAS  CHANGED  PRMT ( 5 )  TO  NON-ZERO. 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
THE  EXTERNAL  SUBROUTINES  FCT ( X * Y * DERY )  AND 

OUTP( X* Y, DERY  * IHLF » NDIM* PRMT)  MUST  BE  FURNISHED  BY  THE  USER. 

METHOD 

EVALUATION  IS  DONE  BY  MEANS  OF  FOURTH  ORDER  RUNGE-KUTTA 
FORMULAE  IN  THE  MODIFICATION  DUE  TO  GILL.  ACCURACY  IS 
TESTED  COMPARING  THE  RESULTS  OF  THE  PROCEDURE  WITH  SINGLE 
AND  DOUBLE  INCREMENT.  1 

SUBROUTINE  RKGS  AUTOMATICALLY  ADJUSTS  THE  INCREMENT  DURING 
THE  WHOLE  COMPUTATION  BY  HALVING  OR  DOUBLING.  IF  MORE  THAN 
10  BISECTIONS  OF  THE  INCREMENT  ARE  NECESSARY  TO  GET 
SATISFACTORY  ACCURACY*  THE  SUBROUTINE  RETURNS  WITH 
ERROR  MESSAGE  3 1  ILF ^ 1 1  INTO  MAIN  PROGRAM. 

TO  GET  FULL  FLEXIBILITY  IN  OUTPUT*  AN  OUTPUT  SUBROUTINE 
MUST  BE  FURNISHED  BY  THE  USER. 

FOR  REFERENCE*  SEE 

RALSTON/WILF,  MATHEMATICAL  METHODS  FOR  DIGITAL  COMPUTERS* 
WILEY*  NEW  YORK/LONDON*  1960  *  PP. 110-120. 
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SUBROUTINE  RKGS < PRMT f Y r BERY f NDIH t IHLF » FCT f OUTP  f AUX > 


C 

C 

C 

C 


C 

c 


c 

c 

c 


c 

c 


c 
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BI MENS I ON  Y  <  1  >  f DER Y < 1 >  f  AUX  <  8  f 1 >  f  A  <  4  >  f  B  <  4  >  f C  <  4  >  f  PRMT ( 1 > 
DO  1  I-l »NDIM 

1  AUX<8fI)=.06A666A7*DERY<I> 

X=PRMT ( 1 ) 

XEND-PRMT  <  2  > 

H=PRMT ( 3 ) 

PRMT  <5)-0.» 

CALL  FCT<  X  f  Y  f DERY ) 

ERROR  TEST 

IF ( H* ( XEND-X ) >38f37f2 

PREPARATIONS  FOR  RUNGE-KUTTA  METHOD 

2  A<1>=.5 

A < 2) =.2928932 
A<3)=1. 707107 
A<4)=.16A66A7 
B  < 1 ) =2 . 

B(2)  =  l . 

B<3)=1 ♦ 

B(4)=2. 

C<1>=.5 
C(2)=. 2928932 
C<3)=1. 707107 
C  <  4 )  = .  5 

PREPARATIONS  OF  FIRST  RUNGE-KUTTA  STEP 
DO  3  1=1 fNDIM 
AUX  < 1 f I ) =Y  < I > 

AUX ( 2 » I >  =DERY  < I > 

AUX(3»I)=0. 

3  AUX ( 6  f I ) =0 . 

IREC=0 

H=H+H 

IHLF=-1 

ISTEP=0 

IEND=0 


START  OF  A  RUNGE-KUTTA  STEP 

4  IF ( (X+H-XEND )#M)7f6f5 

5  H=XEND-X 

6  IEND=1 


7 

8 
9 


RECORDING  OF  INITIAL  VALUES  OF  THIS  STEP 
CALL  OUTP <XfYf DERY f IREC f NDIM f PRMT > 

IF (PRMT (5) )40f8f40 

I  TEST =0  .  . 

ISTEP=ISTEP+1 


rbA' 
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C  START  OF  INNERMOST  RUNGE-KUTTA  LOOP 
J=1 

10  AJ=A< J> 

B  J=B  <  J  > 

CJ=C< J) 

00  11  I-lfNDIM 
R1=H#DERY<I> 

R2=AJ*<R1-BJ*AUX<6»I>> 

Y< I )=Y  < I >+R2 
R2=R2+R2+R2 

11  AUX(6r I )~AUX<6» I )+R2-CJ#Rl 
IF  <  J-4 ) 12 » 15  > 15 

12  J=J+1 

IF ( J-3) 13fl4»13 

13  X=X+.5*H 

14  CALL  FCT(X»YfDERY) 

GOTO  10 

END  OF  INNERMOST  RUNGE-KUTTA  LOOP 


TEST  OF  ACCURACY 

15  IF<ITEST)16»1A»20 

IN  CASE  ITEST =0  THERE  IS  NO  POSSIBILITY  FOR  TESTING  OF  ACCURACY 

16  DO  17  I»1>NDIM 

17  AUX<4»  I )=Y< I ) 

ITEST=1 

ISTEP=ISTEP+ISTEP-2 

18  IHLF=IMLF+1 
X=X-H 
H=.5*H 

DO  19  I»1 »NDIM 
Y ( I ) =AUX ( 1 » I ) 

DERY ( I ) “AUX ( 2  f I > 

19  AUX<&>  I )=AUX(3» I ) 

GOTO  9 


IN  CASE  ITEST-1  TESTING  OF  ACCURACY  IS  POSSIBLE 

20  IM0D=ISTEP/2 

IF( ISTEP-IM0D-IM0D)21 »23»21 

21  CALL  FCT (  X  >  Y » DERY ) 

DO  22  I=1>NDIM 
AUX <5»I)=Y<I) 

22  AUX ( 7 » I >  =DERY  < I ) 

GOTO  9 

COMPUTATION  OF  TEST  VALUE  BELT  * 

23  DELT=0. 

DO  24  1=1 »NDIM 

24  DELT=DELT+AUX ( 8» I ) #ABS<  AUX<  4 » I )-Y < I )  > 

IF (DELT-PRMT  <4> )28»28r25 


A- 7  - 


non  o  o  no 


C  ERROR  IS  TOO  GREAT 

25  IF< IHLF- 10 > 26 » 36 *36 

26  DO  27  1=1 *NDIM 

27  AUX (4*1) =AUX < 5 » I > 

ISTEP=ISTEP+IST£P-4 
X=X-H 
IEND=0 
GOTO  18 

RESULT  VALUES  ARE  GOOD 

28  CALL  FCT <X»Y »DERY> 

DO  29  1  =  1  *NDIM 
AUX  < 1 r I )  -Y ( I  ) 

AUX ( 2 » I ) =  DERY  < I ) 

AUX(3f I )=AUX(6f I ) 

Y ( I ) =AUX ( 5 » I  > 

29  DERY<I)=AUX(7f I) 

CALL  OUTP<X-Hr YrDERY *  IHLF *NDIM*PRMT > 
IF (PRMT < 5) >40*30*40 

30  DO  31  I=1*NDIM 
Y ( I ) =AUX (1*1) 

31  DERY  < I )=AUX<2» I ) 

IREC=IHLF 
IF( I END >32 *32 *39 

INCREMENT  GETS  DOUBLED 

32  IHLF=IHLF-1 
I STEP- I STEP/2 

H 

IF ( IHLF >4*33*33 

33  IM0D=ISTEP/2 
IF  <  I  STEP- 1  MOD*- 1  MOD >4*34*4 

34  IF ( DELT- . 02*PRMT < 4 > ) 35 * 35 * 4 

35  IHLF=IHl.F-l 
ISTEP=ISTEP/2 
H=H+H 
GOTO  4 


RETURNS  TO  CALLING  PROGRAM 

36  IHLF=1 1 

CALL  FCT <X» Y»DERY) 

GOTO  39 

37  IHLF-12 
GOTO  39 

38  IHLF=13 

39  CALL  OUTP<X* Y*DERY » IHLF »NDIM*PRMT> 

40  RETURN 
END 


Vh- 
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